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ABSTRACT 

How do galaxies move relative to one another? While we can examine the motion 
of dark matter subhalos around their hosts in simulations of structure formation, 
determining the orbits of satellites around their parent galaxies from observations 
is impossible except for a small number of nearby cases. In this work we outline 
a novel approach to probing the orbital distributions of infalling satellite galaxies 
using the morphology of tidal debris structures. It has long been understood that the 
destruction of satellites on near-radial orbits tends to lead to the formation of shells 
of debris, while those on less eccentric orbits produce tidal streams. We combine an 
understanding of the scaling relations governing the orbital properties of debris with 
a simple model of how these orbits phase-mix over time to produce a ‘morphology 
metric’ that more rigorously quantifies the conditions required for shells to be apparent 
in debris structures as a function of the satellite’s mass and orbit and the interaction 
time. Using this metric we demonstrate how differences in orbit distributions can alter 
the relative frequency of shells and stream structures observed around galaxies. These 
experiments suggest that more detailed modeling and careful comparisons with current 
and future surveys of low surface brightness features around nearby galaxies should 
be capable of actually constraining orbital distributions and provide new insights into 
our understanding of structure formation. 

Key words: galaxies: kinematics and dynamics - dark matter - galaxies: haloes - 
galaxies: statistics 


1 INTRODUCTION 

In the modern cosmological picture large galaxies are built 
up over time by the hierarchical merging of many smaller 
projenitor systems (e.g. White & Rees 1978). Since the 
galactic luminosity function has a (truncated) power-law 
slope (Schechter 1976, and many others) we expect that the 
vast majority of these mergers will have one parent that 
is significantly more massive than the other. Cosmological 
simulations have shown that Milky Way-mass galaxies have 
typically experienced less than one major merger (with mass 
ratio ^ > 0.1) since redshift z = 1 (Fakhouri, Ma & Boylan- 
Kolchin 2010) and so minor mergers are an important part 
of their total mass accretion rate at late times (Oser et al. 
2010), although possibly subdominant to ‘diffuse’ accretion 
of intergalactic gas (Fakhouri & Ma 2010). 

Observationally, the appearance of debris from minor 
mergers can be broadly divided into two morphological 
categories. Stream-like substructures stretch approximately 
along the progenitor’s orbit, sometimes wrapping around the 
host multiple times, for example as seen in the recent deep 
imaging of NGC 5907 by Martinez-Delgado et al. (2008). 
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Shell-like structures, like those in the vast complex around 
NGC 474 (Due et al. 2013), may extend significantly both 
along and perpendicular to the path of their disrupted par¬ 
ent, forming an umbrella-shaped distribution of stars and/or 
sharp edges in the light distribution, frequently interleaved 
with each other. 

Interest in the detailed study of debris structures has 
been stimulated over the last two decades in part by the 
discovery of numerous streams around the Milky Way, in¬ 
cluding the tidal tails of of Palomar 5, a globular cluster 
(Odenkirchen et al. 2001), the stream of the Sagittarius 
dwarf galaxy (Majewski et al. 2003), the GD-1 stream (Grill- 
mair & Dionatos 2006), and the Orphan stream (Grillmair 
2006; Belokurov et al. 2007). The formation of tidal streams 
is well understood, which makes them powerful tracers of the 
host galaxy’s potential and history: individual stars, heated 
either by tidal forces from the host or internal two-body 
interactions in the case of globular clusters, leave the satel¬ 
lite through the Lagrange points that mark saddle points in 
the the system’s effective potential. Stars that leave through 
the inner Lagrange point have lost energy relative to the 
satellite and have slightly shorter orbital periods, causing 
them to stretch out over time into a leading tail, while par¬ 
ticles that have gained energy have longer orbital periods 
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and lag behind, forming a trailing tail. This model can give 
estimates of the width and length of streams given the or¬ 
bit, progenitor mass and interaction time, and analytic rep¬ 
resentations of this picture have been used to successfully 
produce realistic stream models without full N-body sim¬ 
ulations (Johnston, Sackett & Bullock 2001; Kiipper, Lane 
& Heggie 2012; Sanders 2014; Gibbons, Belokurov & Evans 
2014; Bovy 2014; Amorisco 2015; Fardal, Huang & Weinberg 
2014). 

While there is yet to be entirely convincing evidence 
of shell structures around the Galaxy despite several can¬ 
didates (such as the Triangulum-Andromeda and Hercules- 
Aquila clouds, Rocha-Pinto et al. 2004; Belokurov et al. 
2006; Deason et al. 2013), they are an important class of 
features in external galaxies. Perhaps the original descrip¬ 
tion of substructure as a shell occurs in the Atlas of Pe¬ 
culiar Galaxies (Arp 1966), in reference to Arp 230. The 
Palomar 200-inch plates show at least six shells perpendic¬ 
ular to the principal axis. Many other galaxies in the At¬ 
las have shell-like debris structures, e.g. Arp 223 and 227, 
an early hint that they are common - although some pro¬ 
cessing may be necessary to reveal them due to their low 
surface brightness (Malin & Garter 1983). Simulations of 
shells followed, starting with restricted n-body models where 
the host potential was static and the disrupting satellite’s 
self-gravity was assumed negligible (Quinn 1984; Dupraz & 
Gombes 1986). Subsequent work modeled the satellite po¬ 
tential self-consistently (Piran & Villumsen 1987; Heisler 
& White 1990). Both showed that radial mergers can re¬ 
produce the observed shell properties. Observational results 
disfavor models where the shells are formed through the 
production of density waves in the stellar outskirts of the 
host galaxy by a perturber (the Weak Interaction Model, 
Thomson & Wright 1990) as an alternative to the merger 
scenario in many systems (e.g. Turnbull, Bridges & Garter 
1999; Wilkinson et al. 2000; Schiminovich, van Gorkom & 
van der Hulst 2013). More recent theoretical treatments of 
shell formation have focused on recovering the host gravi¬ 
tational potential from the line-of-sight velocity distribution 
of the debris (Merrifield & Kuijken 1998; Ebrova et al. 2012; 
Sanderson & Helmi 2013). Detailed models can capture the 
shell density distribution for purely radial orbits (Sanderson 
& Bertschinger 2010) but no simple description of the condi¬ 
tions under which minor merger debris forms a shell rather 
than a stream has emerged until very recently (Amorisco 
2015, this work). 

Overall, previous studies demonstrate that merger de¬ 
bris’ properties are functions of the interaction time (length 
of streams or angular extent and number of shells) and mass 
ratio (which influences the spread of the debris, both in posi¬ 
tion space and energetically). Since the disruption process is 
well enough understood to actually be invertible (e.g. John¬ 
ston 1998; Helmi & White 1999), the population of streams 
with different extents and surface brightnesses can conceptu¬ 
ally be used to ascertain the rate of minor mergers of differ¬ 
ent mass ratios (as outlined in Johnston, Sackett & Bullock 
2001; Johnston et al. 2008) and place more stringent con¬ 
straints on structure formation than a coarser classification 
as major or minor mergers. However, the signatures of most 
accretion events are expected to be extremely low surface 
brightness (LSB), starting near 28 mag arcsec”^, with the 
majority below ~ 30 mag arcsec”^ (Johnston et al. 2008; 


Gooper et al. 2010); fortunately, exploration of the LSB 
universe is rapidly advancing with the availability of wide- 
held cameras, specialized instruments, and LSB-optimized 
observing techniques (Janowiecki et al. 2010; Miskolczi, 
Bomans & Dettmar 2011; Atkinson, Abraham & Fergu¬ 
son 2013; Martinez-Delgado et al. 2010; Abraham & van 
Dokkum 2014; Due et al. 2015). The future Large Synoptic 
Survey Telescope (LSST) is expected to image the Southern 
sky with surface brightness sensitivity 5 magnitudes deeper 
than current wide-held surveys such as the Sloan Digital Sky 
Survey (SDSS) (Ivezic et al. 2008). 

The promise of these current and near-future investiga¬ 
tions motivates a reexamination of the utility of substruc¬ 
ture as probes of galaxies’ pasts. Since the debris structures 
that result from accretion events with near-radial infall are 
morphologically distinct from those produced along less ec¬ 
centric orbits, substructure counting could conceptually be 
sensitive to the accretion rates as a function of mass, time 
and orbit - adding a new dimension of information about 
the accretion history that is otherwise difficult to access. 
It is not possible to reconstruct the full orbits of most ob¬ 
jects from observations since only the projected distance and 
line-of-sight component of motion can be measured outside 
the Local Group and the typically small number of known 
galaxy-mass satellites in individual systems (outside clus¬ 
ters) precludes the use of statistical tools such as the Jeans 
equations without stacking (Herbert-Fort et al. 2008). While 
each galaxy may have only one or fewer debris structures ap¬ 
parent, a well-constructed survey of galaxies to faint surface 
brightness could provide constraints on the collective dis¬ 
tribution of infalling objects. A strong grasp on the orbital 
distribution of accreting objects provides an important next 
step towards reducing the number of degrees of freedom in 
the choice of cosmological models, similar to the way that 
the merger fraction as measured by galaxy pairs has been 
used to restrict combinations of the principle cosmological 
parameters (Garlberg 1991; Gonselice et al. 2014), and might 
have implications for our understanding of how baryons in¬ 
habit low-mass haloes. 

In this paper we explore the idea that debris structures 
carry useful information about the satellite orbital infall dis¬ 
tribution. In Section 2, we briefly describe N-body simula¬ 
tions used to examine debris formation throughout the rest 
of this work. In Section 3 we extend the simple stream¬ 
building picture described above to also encompass shells 
and use the results to define the morphology metric /i. In 
Section 4 we use our metric to demonstrate that the pop¬ 
ulation of shells and streams is indeed sensitive to the in¬ 
falling subhalos’ orbital velocity distribution and that these 
populations could potentially be used to provide interesting 
constraints on cosmological models that predict different or¬ 
bital infall distributions, presuming a variety of intermediate 
challenges can be overcome. We consider the observational 
and modeling advances necessary to reasonably apply this 
technique in Section 5. Section 6 concludes. 


2 N-BODY SIMULATIONS 

We performed a set of dark-matter-only (but see 4.2, where 
baryonic effects are considered) N-body simulations with the 
self-consistent held basis function expansion code (Hernquist 
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& Ostriker 1992) to explore the formation of debris struc¬ 
tures across a wide range of orbital and galactic parameters. 
In each simulation, a 10® particle NFW-profile (Navarro, 
Frenk & White 1996) satellite was inserted at the apogalac- 
ticon of its orbit in a static, NFW host potential. The satel¬ 
lite evolves first in isolation to ensure that it equilibrates 
and then the host potential is gradually turned on over 10 
satellite internal dynamical times to reduce artificial gravi¬ 
tational shocking. Total energy is conserved to ~1% of the 
satellite internal potential energy during the 8 Gyr integra¬ 
tions. The simulation is not halted if the satellite is com¬ 
pletely disrupted since the debris continues to evolve. 

The satellites are initialized with masses m/M© = 6.5 x 
10®, 6.5 X 10^, 6.5 X 10®, and 6.5 x 10®, where m is the mass 
enclosed at 35 NFW scale radii, out to which particles were 
realized. The scale radius ro was adjusted for each mass so 
that the mean density inside the scale radius is the same for 
each, with a value of ro = 0.86 kpc for the 6.5 x 10® Mq 
satellite. This scaling gives them the same fractional mass 
loss rate along a given orbit. 

We use a host that is broadly consistent with expecta¬ 
tions for a Milky Way-scale dark matter halo, choosing a 
viral mass of 1.77 x 10^® M©, virial radius of 389 kpc and a 
scale radius of 24.6 kpc. We chose a set of orbits with total 
energy equal to that of circular orbits at 25, 50, and 100 
kpc in this potential and varied the angular momentum in 
twelve steps between 0.05 and 1.0 times the angular momen¬ 
tum of the circular orbit for each energy, denoted by Tcirc- 
The lowest angular momentum case for these orbits is quite 
radial; the rcirc = 25 kpc, L/Lcirc = 0.05 simulation has 
a perigalacticon distance of ~ 0.65 kpc and an apogalacit- 
con of ~ 40 kpc. This orbit is a reasonable match to the 
extent of the shell systems seen around NGC 4651 (Foster 
et al. 2014) and MGG-5-7-1 (Schiminovich, van Gorkom & 
van der Hulst 2013), both of which have dark matter halo 
masses estimated to be of the same order as our selected 
value. 

These simulations were used to check the scaling re¬ 
lations in Section 3.2 and to substantiate a choice of the 
morphology metric that will separate shells from streams. 


3 RESULTS I: DEFINING A MORPHOLOGY 

METRIC 

To build a statistic that will allow us to quickly assess what 
debris from a disrupted satellite will look like given a par¬ 
ticular set of merger parameters we first presume that the 
host potential is sufficiently well known (at the population- 
statistics level) that the fact that a shell forms at all places 
a constraint on the characteristics of an individual merger. 
This is in contrast to previous work (Merrifield & Kuijken 
1998; Ebrova et al. 2012; Sanderson & Helmi 2013) that 
uses shells to constrain host halo parameters. Those analyses 
concentrated on modeling shell systems around individual 
galaxies in detail through density and line-of-sight velocity 
matching, while here we are interested in the general con¬ 
ditions that lead to shells and streams (without assuming 
any spectral information) in order to examine the implica¬ 
tions of observed merger morphologies for the distribution of 
satellite orbital properties around a large sample of galaxies. 

In this section we begin by reviewing the properties of 
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Figure 1. A visualization of the angles a (black dashed lines) 
and Al/) (gray dashed lines) for orbits of varying eccentricity. The 
solid black lines in each panel show an orbit with energy equal 
to the 25 kpc circular orbit but with the specified fraction of 
that orbit’s angular momentum. The dotted circle marks the half¬ 
period radius; a particle will spend equal time inside and outside 
this radius. For the orbits shown it varies from 25.7 kpc for the 
most circular orbit to 31 kpc for the most eccentric orbit. Inside 
a a particle is closer in time to apogalacticon than perigalacticon. 


test particle orbits in extended mass distributions, then use 
these insights to understand what properties of mergers lead 
to the creation of shells and streams (Section 3.1). Identify¬ 
ing the debris’ energy and angular momentum dispersions 
(Te and ul as the quantities of interest in the morphological 
distinction, we describe a method to estimate the disper¬ 
sions in Section 3.2. Finally we fold this understanding into 
a morphology metric that can predict the debris morphology 
in Section 3.3. 


3.1 Properties of orbits 

Knowledge of the orbits of test particles in simple poten¬ 
tials is an important first step towards understanding merg¬ 
ers in their entirety, as shown in the original classic work 
by Toomre & Toomre (1972). In the case of debris result¬ 
ing from a dwarf galaxy entering a Milky Way-size galaxy’s 
halo we can begin by assuming that the remnant’s grav¬ 
itational infiuence on the unbound particles is negligible. 
This assumption is well justified for mass ratios smaller than 
^ = m/Mhost = 10“"^ (Ghoi, Weinberg & Katz 2007) and in 
our simulations the debris’ conserved quantities are unaf¬ 
fected by the satellite after unbinding until ^ > 10“®. Addi¬ 
tionally, we expect the host potential at the typically large 
radii considered below to be dominated by the host’s (spher¬ 
ical) monopole term, since any effect from disks must be 
subdominant given that both streams and shells are found 
in galaxies with elliptical morphologies as well. While devi¬ 
ations from spherical symmetry will introduce a variety of 
interesting and potentially important effects such as orbital 
plane precession and chaos, as a limiting case spherical hosts 
are useful for building intuition. 

Orbits in spherical potentials can be uniquely identi¬ 
fied by their conserved quantities - the total energy E and 
total angular momentum L - up to the orbital plane orien¬ 
tation, which we take as the x — y plane for convenience. 
These, combined with the potential parameters, determine 
the shape and properties of the orbit: the radii of the turn¬ 
ing points at apogalacticon and perigalacticon, radial orbital 
period Tr, and the precession per orbit between apogalactica 
A-)/). Fig. 1 illustrates how orbits change when the circular¬ 
ity L/Lcirc is varied. The precession angle and perigalac- 
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Figure 2. Various orbital descriptors as a function of energy and angular momentum for test particles orbiting in the spherical NFW 
host halo used for the N-body simulations (Section 2). White space indicates regions that are inaccessible in this host potential. The 
periods and apogalacticon distance are strong functions of energy but weak functions of angular momentum, whereas the changes in 
precession and perigalacticon distance are dominated by angular momentum. The “width” a is primarily determined by L/Lcirc- 


tic distance decrease with decreasing circularity while the 
apogalactic distance increases. We will also find it useful to 
consider the angular ‘width’ of a single petal of the rosette 
traced by an orbit, which varies systematically with angu¬ 
lar momentum; near-radial orbits spend the vast majority 
of their time near apogalacticon, covering only a small an¬ 
gle in azimuth during a large fraction of their radial period, 
while near-circular orbits may cover hundreds of degrees in 
azimuth during the same portion of radial phase. To charac¬ 
terize this difference, we define a as the angle through which 
a particle moves during the outer half of its radial period. 
This is a proxy for the angle subtended by the orbit’s minor 
axis, which does not exist for these non-closed orbits. 

Fig. 2 illustrates how the orbital properties vary as 
a function of energy and angular momentum. Each panel 
shows the available space for bound orbits with energies less 
than that of a circular orbit at 150 kpc in the same host 
as described in Section 2 - an NFW halo with virial mass 
1.77 X 10^^ Mq and a scale radius of 24.6 kpc. The radial 
period Tr is a strong function of energy but only weakly 
dependent on angular momentum, while for the precession 
per orbit the opposite is true. Apogalacticon and peri¬ 
galacticon distances are also primarily functions of energy 
and angular momentum, respectively, but the dependencies 
are weaker here than for Tr and A-i/i. The angle a follows 
neither of these trends but instead depends on the orbit’s 
circularity rather than the absolute values of E or L. 

Fig. 3 shows the visible effects of sampling the prop¬ 
erties shown in Fig. 2 across a small range in the {E, L) 
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Figure 3. Variation of test-particle orbits due to small changes 
in energy (left) or angular momentum (right); points show the 
positions of particles that are displaced from from the gray central 
orbit in each quantity by ±2.5% (±10%). Both figures use the 
same reference orbit and are integrated in the same potential 
and points are displayed at the same times. While variation in 
E displaces particles primarily along the reference orbit, varying 
L causes precession that distributes particles in azimuth with a 
range in position angle that grows with time. 


parameter space. The gray rosette represents the position of 
a reference orbit at every time step of a simple numerical 
leapfrog integration in the standard halo described above; 
it has energy equal to that of a circular orbit at 37 kpc but 
only 30% as much angular momentum, i.e. L/Ldrc = 0.3. 
The points near the rosette show the positions of a few or¬ 
bits with a distribution of properties around the reference 
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Figure 4. Debris formation resulting from simple test particle models for the satellite (left) and and N-body simulations for the same 
(right). The background potential has a static NFW profile. Individual particles are color coded by either their energy (top row) or 
angular momentum (bottom row) to illustrate how offsets in orbital quantities dictate their final spatial position. 


orbit at a few different times. Ten particles start at (x,y,z) 
= (25,0,0) kpc and for the left (right) panel their energies 
(angular momenta) were varied by ±2.5% (±10%) in equal 
spacings. The energies and angular momenta set the magni¬ 
tude and direction of the initial particle velocities and their 
positions are evolved with the same leapfrog scheme. As ex¬ 
pected from Fig. 2, varying the orbital energy while holding 
L constant primarily affects the orbital periods of particles, 
spreading them out in radial phase but keeping their or¬ 
bits otherwise similar for cosmologically significant interac¬ 
tion times (here 4 Gyr), giving them the appearance of a 
stream tracing the ‘progenitor’ orbit. The most obvious de¬ 
viation from the reference orbit is an increased apogalacticon 
distance for the trailing particles and decreased apogalacti¬ 
con distance for the leading ones; this also results in the 
rosette petal tips, as traced by the debris, becoming more 
rounded than the reference orbit as high-energy particles 
move farther out when they are near apogalacticon. In con¬ 
trast, when angular momentum is varied but energy is held 
constant particles remain tightly clustered in radial phase 
because their periods are nearly equal. However, they slowly 
spread azimuthally due to differences in the rate at which 
their apogalactica precess, resulting in a ‘shell’ of particles 
at the same radius near apocenter since the apogalactic dis¬ 
tance varies only weakly with L. The selected values of the 
conserved quantities and halo do not change the qualitative 
effect. 

In reality, streams and shells will contain particles offset 
from the satellite in both energy and angular momentum. 
The left panels of Fig. 4 show the result of integrating 10® 
test particles on the same orbit as Fig. 3 but where the parti¬ 
cles are given independent Gaussian distributions in energy 
and angular momentum, with dispersions cte and ctl respec¬ 
tively. The energy dispersion is 2% of the mean orbital en¬ 


ergy in all panels. The particles are color-coded by their 
energy (top panels) or angular momentum (bottom panels), 
with blue indicating the lowest values and red the highest. 
This simple setup produces quite reasonable ‘debris’ (com¬ 
pare with the full N-body simulations in the right panels of 
Fig. 4), which is perhaps surprising, but this model is es¬ 
sentially a minimum-complexity version of the phase-space- 
distribution method to produce streams (Kiipper, Lane & 
Heggie 2012; Gibbons, Belokurov & Evans 2014; Amorisco 
2015). As before the qualitative result of energy and angu¬ 
lar momentum sorting is independent of the specific orbit or 
values of ite and ctl chosen. 

Fig. 4 also shows that when the fractional variation in 
E and L is equal, (ctl/L)/(cte/E) = 1 (left panels), the dis¬ 
tribution has a clearly stream-like appearance. Particles of 
the same energy (i.e. same color in upper left panel) have 
very nearly the same period, reaching apogalacticon at the 
same time. The precession rate, dominated by variation in 
angular momentum, spreads the particles in azimuth and 
creates the sharp, rainbow-colored edge at the ends of the 
rosette petals (lower left panel). The precession is insuffi¬ 
cient to significantly alter the appearance that the debris 
follows the primary orbit. When ((Jl/L)/((Te/E) = 4, how¬ 
ever, the precession dominates and the debris appears as a 
shell because the angle subtended by the radial edge due to 
differences in the precession rate exceeds the width of the 
rosette petals by a factor of ~ 3. 

From previous work it is well known that shells are more 
frequent in the aftermath of radial mergers. Here we have 
additionally demonstrated that, for a given orbit, the ratio 
of fractional dispersions in E and L inside a debris struc¬ 
ture is an important factor in determining the morphology. 
However, cte and ctl are not free parameters but are instead 
set by the orbit and the internal properties of the satellite. 
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As shown in Section 3.2 below and Fig. 5, the actual or¬ 
bital parameter distribution in the debris a) is bimodal, b) 
has a specific, quite non-Gaussian shape in each mode, and 
c) has significant covariance especially for low eccentricity 
orbits (see also their projections into action-angle and fre¬ 
quency space in Sanders & Binney 2013; Bovy 2014; Fardal, 
Huang & Weinberg 2014). This covariance is why the left 
panels of the test panel integration in Fig. 4 have a more 
pronounced apogalacticon edge than real streams: the E- 
L covariance will only allow a small subset of angular mo¬ 
menta at a given energy, and therefore smaller differences 
in integrated precession angle at constant E. We find that 
accounting for this covariance is unnecessary to predict the 
morphology from N-body simulations (Section 3.3) 


3.2 Debris scaling relations 


The typical scales in debris (which characterize the disper¬ 
sions (Te and ctl) can be computed following previous work 
(e.g. Johnston 1998; Binney & Tremaine 2008). Consider a 
satellite moving with velocity Vp at its perigalacticon, where 
its distance from the host galaxy is Rp. Most mass loss will 
occur near pericenter so we anticipate that the debris’ dis¬ 
persions in E and L are set by the relevant satellite scales 
there. Particles are very likely to leave the satellite through 
the effective Lagrange points where the effective potential 
has a saddle point; these points are in general not equally 
spaced but are colinear with the host and satellite centers. 
We approximate their location relative to the satellite cen¬ 
ter by the tidal radius rude where the gravitational forces 
would balance if the galaxies were point masses on circular 
orbits: 

/ ^ \ 1/3 

= [sMiRp)) ^ 


where m and M{Rp) are the mass of the satellite and the 
mass of the host enclosed inside Rp, respectively. 

We take energy scale Ca of the debris to be the differ¬ 
ence in the host’s gravitational potential energy across the 
satellite diameter at pericenter using the simple linear ap¬ 
proximation 


Cs — 2rtide 


a-F 

M 


Rp 


( 2 ) 


while the angular momentum scale is computed from 


la = AL = AVR + VAR. (3) 


Using the satellite internal velocity dispersion a = 
^/Gm/r\uJ for AV, AR = 2rtide, and Vp = ^JGM/Rp, 
one finds 


la — crRp -\- 2Vprtide — (V^ -t- 2)sL (4) 


Thus Ba and la are specified by the orbital and galactic pa¬ 
rameters and therefore, given our conclusions in Section 3.1, 
so is the debris morphology; varying Ca and la independently 
for a given satellite and orbit as in the test particle integra¬ 
tion of Fig. 4 is unphysical. 

The efficacy of these scales is demonstrated in Fig. 5. 
Shown are the results of the simulations described in Sec¬ 
tion 2 for Ccirc = 25 kpc with the m = 6.5 x lO^M© run in 
red and m = 6.5 x IO^Mq in blue, after 8 radial periods (ap¬ 
proximately 4.3 Gyr). In each panel satellite particles are 



Figure 5. Energy and angular momentum offsets of bound and 
unbound particles. 


plotted in {E, L) space as offsets from the satellite center 
of mass values Eorb and Lorb', all but the bottom-left panel 
have their axes rescaled by Cs and L- In the upper panels 
black points indicate particles that are still bound to the 
satellite while colored points have become unbound through 
tidal stripping. The color tone (lightness) of colored points 
represents the angular momentum of the satellite orbit they 
come from, with darker tones indicating more circular orbits. 
Snapshots of the disruption near apocenter (top row) and 
pericenter (middle row) for a range of angular momenta are 
provided; the bound particle distribution rotates about the 
satellite center {Eorb,Lorb) between these directions, filling 
the space with unbound particles (outside an exclusion zone 
defined by the satellite’s interaction with its own debris). 
The energy and angular momentum of individual particles 
is constant to high precision after they have been unbound 
for a short time, consistent with the assumption of no self¬ 
interaction at these mass ratios (^ ~ 10“® — 10“®) 

In the bottom row, the unbound particles from sim¬ 
ulations with a wide range of satellite angular momenta 
L/Lcirc = 0.1, 0.2 ... 0.9 are plotted together for the two dif¬ 
ferent mass satellites. Shown are absolute offsets (left) and 
rescaled offsets (right). Although there is some remaining 
angular momentum dependence the scalings do an excellent 
job of removing the mass dependence and reducing the vari¬ 
ation in median offsets across circularities. While corrections 
could be made to ea and L to account for the remaining an¬ 
gular momentum dependence and other effects such as mass 
loss over time, dynamical friction, the changing in tidal ra¬ 
dius as a function of instantaneous galactocentric distance, 
or more general potentials, we consider these simple esti¬ 
mates sufficient for our purpose. 
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Figure 6. Time evolution of debris. Left: the simulation panels show snapshots from the rcirc = 25 kpc, L/Lcirc = 0.7 simulations of 
satellites with m = 2.5 X IO^Mq (top row) and m = 2.5 X 10*Mq (bottom row) at selected times, viewed along the axis perpendicular 
to the orbital plane.. Right: the various angles, corresponding to the displayed simulations, that contribute to the morphology metric. 


3.3 A morphology metric 


So far we have established that (i) modest variations in or¬ 
bital parameter space have effects on properties (orbital pe¬ 
riods, precession) that are important to the shape of debris 
on the timescale of satellite disruptions, (ii) which proper¬ 
ties will be affected by such variations depend strongly on 
whether it is E or L that is varied, and (hi) satellite dis¬ 
ruption produces debris in a predictable section of orbital 
parameter space in the vicinity of the satellite’s own orbit. 

Using these facts, we propose a metric for estimating 
whether a debris structure would be visually classified as a 
shell or a stream by comparing the contribution of energy- 
dominated and angular momentum-dominated effects to its 
shape. As illustrated in Fig. 4, a reasonable choice for the 
factor determining morphological assignment is whether or 
not the shell edge subtends a larger range of position angle 
than a single rosette petal. We therefore need a measure of 
1) the rosette petal angle, already defined as a, 2) the angle 
through which debris has spread due to differences in az¬ 
imuthal precession rate as a function of time, l and 3) the 
angle subtended due to stream-like orbital period variations 
in the debris, IHb, if less than a whole petal. 

The near-complete separability of effects from E and 
L allows a simple linear approximation. Given a satellite- 
host pair and an orbit {Eorb, Lorb), the angular size of the 
constant-energy edge caused by precession of apogalactica 
is just the differential precession per orbit with respect to L 
at Eorb times the angular momentum scale and the number 
of orbits Aorb: 


\&L = IsNoih 


dAip 

~dr 



( 5 ) 


where the contribution proportional to dA'tp/dE is neglected 
due to the weak dependence of A'i/) on E. Similarly, since the 


radial period is a weak function of L the debris’ extent from 
energy variations can be approximated by considering the 
added (subtracted) angle that short- (long-) period, low- 
(high-) energy debris moves through, given by 


'i’E = CsNoih 


Alp dTr 

~T7~dE 



( 6 ) 


We find that e > l generically; dynamically young de¬ 
bris will always be stream-like. We restrict to be less 
than a so that we are considering individual rosette petals 
(see Fig. 1) and so we define our morphology metric /r as 


where 


r- = 




( 7 ) 


\1/b = min (a, tks) . (8) 

Fig. 6 illustrates the time evolution of the angles tk e and 
'i’L StS computed from Equations (5) and (6) for two mergers 
in our standard host halo with rcirc = 25 kpc and circularity 
L/Lcirc = 0.7 but differing masses. When /r is larger than 
1 shells should become apparent in the debris. The rosette 
petal width a (horizontal dotted line) is a function of orbit 
only and therefore constant with respect to ^ and time. The 
stream angle (solid black line) is larger than ^kn (dashed 
black line) at early times; both angles increase linearly so 
fi{t) is constant until \k_B exceeds a. Eventually \k_B > a and 
so tklj = a thereafter. Since the integrated precession angle 
\k L continues to grow, the morphology becomes increasingly 
dominated by angular momentum effects until Vkn = Vk);, i.e. 
/r = 1 (vertical black dotted line), when the merger has pro¬ 
duced a shell. At the last snapshot both mergers show shells 
as expected from to the moderately eccentric orbit but the 
mass dependence results in a significantly faster transition 
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Figure 7. Left: the theoretically motivated divider = 1 (black line) effectively separates debris structures that have been visually 
classified as containing shells (red circles) from those marked as streams (blue triangles). Most objects that could not be confidently 
assigned were dynamically young. Some objects with > 100° have been omitted for clarity. Right: selected objects with fi = 10.6, 
2.8, 1.4, 0.49, 0.45, and 0.33, viewed along the axis perpendicular to the orbital plane. 


for the higher mass case, by more than a factor of two in this 
example. The time dependence of /r adequately captures the 
morphological evolution. 

In this model all debris structures inevitably evolve into 
shells. At late times phase mixing will reduce the density en¬ 
hancements to the point where debris may not be identified 
as such and it might be desirable to exclude mergers that 
have progressed to this stage from our accounting in order 
to fairly compare with observations; however, this process is 
closely tied to related issues of surface brightness and orien¬ 
tation effects that are left for future work. 

To further test the morphology metric, 200 snapshots 
were selected at random from the 14,400 outputs of the sim¬ 
ulations described in Section 2 and were morphologically 
classified by eye, viewed along the axis perpendicular to the 
orbital plane. Equations (5) and (8) were then used to com¬ 
pute the energy- and angular momentum-induced position 
angle variations from the simulation initial conditions. Fig. 7 
shows that debris marked as containing shells versus those 
that appear as only streams are cleanly separated in (v& n, 
space, and that /r = 1 can effectively divide the differ¬ 
ent structures. In the right-hand panels, six representative 
snapshots are shown with a wide range of flfL, and /r, 
illustrating how sorting by decreasing values of the mor¬ 
phology metric shifts from classical shell systems (A, B) to 
“umbrellas” (C) and then streams (D, E, F). 


4 RESULTS II: THE INFLUENCE OF 
ORBITAL DISTRIBUTIONS ON THE 
FREQUENCY OF SHELLS 

The rapid acceleration of computational speed, storage, and 
techniques now allows cosmological N-body simulations of 
structure formation to evolve representative volumes of the 


Universe to 2 ; = 0 in order to track statistical samples of 
haloes down to dwarf galaxy scales, M 200 ~ IO^^Mq (for ex¬ 
ample, Millennium-II, Boylan-Kolchin et al. 2009; Bolshoi, 
Klypin, Trujillo-Gomez & Primack 2011). Of particular in¬ 
terest for our study is the haloes’ accretion history. Typically 
the accretion history is projected into two spaces: the mean 
merger rate dNin/d^/dz(5, M, z) (which has units of mergers 
per halo, for a halo mass M, per unit redshift z per unit mass 
ratio and the orbital inf all distribution which is usually 
described by the probability distribution P(Vr, Ve | M, z) 
for finding satellites with radial and tangential velocities in 
the range Vr + dVr and Ve -I- dVe. 

The morphology metric /r = z, E, z)) 

described in Section 3.3 allows us to make a direct con¬ 
nection between the orbital infall distribution and the ob¬ 
served numbers of shells and streams. The expected number 
of mergers that produce debris structures exhibiting shells 
around a galaxy of mass M can be computed, after relating 
the orbital parameters to their equivalent velocities through 
the assumed halo potential {E,L) !—>■ ^{yr,Ve), by integrat¬ 
ing over the distributions 


-*max / smax 


Yahell(M) = 






dNm 

dCdz"^ (9) 


P{Vr,Ve I ^,M,z)}i{fi-^it)dVrdVed^dz. 


Here I4sc is the host escape velocity, /it is the value of the 
morphology metric chosen to demarcate the transition from 
streams to shells, and H(a;) is the Heaviside step function, 
equal to 1 if a; > 0 and 0 otherwise. Computing Aatream(A/) 
instead simply requires negating the argument of H. The 
quantities NsheiJM) and Aistream(Af) are the model pre¬ 
dictions for the absolute number and relative frequency of 
mergers that create each type of debris structure. It is im¬ 
portant to note that Aaheii counts the number of events that 
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produce shells, not the number of individual surface bright¬ 
ness edges resulting from a merger, and therefore there is 
risk of observational confusion when multiple progenitors 
produce debris around the same host galaxy. Additionally, 
it is assumed that all events inside the integration limits 
are detectable without regard to uncertainties introduced 
by the host light distribution and any survey sensitivity ef¬ 
fects. Nevertheless, this analytic approach allows us to make 
some first estimates of how sensitive the population of shells 
and streams might be to the orbital infall distribution. Sec¬ 
tion 4.1 outlines the current state of measurements of vari¬ 
ous properties of haloes in cosmological simulations. In Sec¬ 
tion 4.2 we implement an approximate representation for 
the luminous matter in the satellites and effective surface 
brightness limits. Combining the two. Section 4.3 makes pre¬ 
dictions for the observed debris population. 



4.1 Parameter and model choices 

We presume that at some future date a merger rate, NFW 
halo concentration relation, and orbital infall distribution 
will all be consistently measured from a single cosmological 
simulation. To schematically illustrate how these predictions 
could be compared to observations we will use some com¬ 
monly cited fits from separate simulations, each of which 
necessarily has its own cosmological parameters, cosmic vari¬ 
ance, systematic errors and resolution effects to consider. 
Any other choice for each may be trivially substituted. 


4-1.1 Merger rate 

The per-halo merger rate provides the overall scaling of the 
absolute number of substructures. Fakhouri, Ma & Boylan- 
Kolchin (2010) combined the wide Millennium-I (Springel 
et al. 2005) simulation with the better mass resolution of 
Millennium-II (Boylan-Kolchin et al. 2009) to achieve good 
statistics on mergers over 5 decades of host mass and merger 
mass ratio from the present to z « 15. They provide a fitting 
function to the mean merger rate per halo dA^m/d^/dz as a 
separable function of host mass M, mass ratio ^ and redshift 


dAr,„ 

d^dz 




f M 


VlOi^Mo 


4 exp 



(l + z)« 


( 10 ) 

where = (0.133,-1.995,0.263,0.0993), (A,^) = 

(0.0104,9.72 X 10“^). The rate per halo is a weak function of 
host mass and redshift but a much stronger function of mass 
ratio; the quality of satellite and host masses measured in 
the observations will have a large effect on the uncertainty 
in When- 


4-1-2 Host halo parameters 

As indicated by Equation 2 and also shown by Amorisco 
(2015), the slope of the host halo’s potential at perigalacti- 
con directly affects the spread in energy and hence the rate 
of shell formation. For spherical NFW haloes of a given mass 
this slope is determined by the concentration parameter 
Cvir = rvir/rs, where rvir is the host virial radius and rg is the 
NFW scale radius. Navarro, Frenk & White (1996) showed 
that there is a tight correlation between concentration and 



Figure 8. Outside-in stripping of mass. Top: “baryons” repre¬ 
sented by the particles with initial binding energies in the 90th 
percentile and higher, are protected from tides until ~ 80% of the 
total mass has been lost independent of the mass, orbital energy, 
and orbital angular momentum of the satellite. Center: Shells in 
dark matter (black) vs baryons (red). Bottom: The baryons (red) 
have smaller energy and angular momentum offsets because they 
are only stripped when the satellite has lost a large fraction of its 
mass (compare the bottom-left panel of Fig. 5) 
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mass; recently Dutton & Maccio (2014) found the following 
relation holds for relaxed haloes with 10^® < M/Mq < 10^® 
to redshift z = 5 with a simulation using the cosmological 
parameters as measured by the Planck satellite (Planck Col¬ 
laboration et al. 2014): 

logic Cvir = a + b logio(M/10^^/i“^M©) 

a = 0.537 -b (1.025 - 0.537)exp(-0.718«^ °®) (U) 

6 = -0.097-b 0.024«. 

4.1.3 Orbital parameter distributions 

The distribution of orbital parameters may depend on a 
number of merger properties: the host mass, the mass ra¬ 
tio, and the redshift of infall are all likely candidates besides 
the chosen cosmology. Past work has disagreed about which 
of these are important. Some suggest ^ dependence (Tormen 
1997; Vitvitska et al. 2002; Jiang et al. 2015) while others 
do not (Gill et al. 2004; Khochfar & Burkert 2006; Wetzel 
2011). Some find redshift dependence (Benson 2005; Wet¬ 
zel 2011) while others who investigated exclude it (Zentner 
et al. 2005). Whether the host mass has any influence is 
also disputed, however there are large variations in dynamic 
range between these studies which may obfuscate this. The 
authors who investigated correlation between orbital param¬ 
eters typically found it (Tormen 1997; Benson 2005; Wang 
et al. 2005; Jiang et al. 2015). 

We consider possible correlation between Vr and Ve 
an important factor since it is required e.g. to capture the 
situation where infalls have a preferred energy. With this 
in mind we will use the findings of Benson (2005) and 
Jiang et al. (2015) to test whether the number of shells and 
streams can distinguish between infall models and there¬ 
fore, given sufficient grips on what range of M, and z 
one is probing, the consistency of simulation and obser¬ 
vation. Benson (2005) provides fits to the probability dis¬ 
tribution P(Vr,Ve I z) for zG {0.0,0.5,1.0} while the re¬ 
sults of Jiang et al. (2015) can be transformed to com¬ 
pute P(V,,Ve I M,C) with M/M© e {10l^l0l^l0l^} and 
5 G {[0.0001,0.005], [0.005,0.05], [0.05,0.5]}. The shapes of 
these infall distributions in {Vr,Ve) space are shown in the 
left panels of Fig. 11 and Fig. 12 

4.2 Converting dark matter simulations to LSB 
features 

4.2.1 Two-component satellites 

Thus far we have considered the disruption of only the satel¬ 
lites’ dark component. The baryonic part is situated deep in 
the satellite potential and is therefore significantly more re¬ 
sistant to tides than the more extended dark matter. This 
can lead to substantial differences in the spatial distribu¬ 
tion of the two components simply because the satellite will 
have lost much of its mass before any stars are unbound 
and therefore the area covered by stars in [E, L) space will 
smaller than computing Cs and P from the initial mass would 
suggest. We represent the fraction of the binding energy dis¬ 
tribution that contains stars by fb. The true value of this 
fraction is not well understood and depends on the details of 
the galaxy formation process but previous work has shown 
that tagging the 10% most bound dark matter particles 


with stars reproduces many properties of the stellar halo 
in high-resolution resimulations of Milky Way-sized galaxies 
(De Lucia & Helmi 2008) and captures the observed surface 
brightness, velocity dispersion, and luminosity distributions 
of Local Group dwarfs modeled from semi-analytic initial 
conditions (Bullock & Johnston 2005) although a smaller 
fraction (1-3%) may be required to ensure that the size- 
luminosity relation of Local Group dwarfs is also respected 
when the halo is live (Gooper et al. 2010). 

Fig. 8 explores the most extreme two-component sce¬ 
nario, where all particles with initial binding energy in the 
90th percentile and above are considered ‘baryons,’ fb = 0.1. 
In the top panel the fraction of baryon-labeled particles that 
are still bound is plotted as a function of the total bound 
mass fraction at each output for 72 of our simulations, those 
with L/Lcirc <0.5; this cut is used to ensure sufficient num¬ 
ber statistics. As expected the baryons are well isolated from 
tidal forces and virtually none of them are removed until 
~ 80% of the initial mass is stripped. This result matches 
that found previously by Villalobos et al. (2012) and Ghang, 
Maccio & Kang (2013). The correlation between the remain¬ 
ing fractions of the two components is very similar across all 
simulations and the deviations that exist do not appear to 
be related to mass ratio, circularity or orbital energy. 

This multiple component model has two important ef¬ 
fects on the morphology metric. First, there is a time de¬ 
lay introduced by the requisite halo stripping; this can be 
less than an orbital period if the orbit is highly eccentric 
or more than a Hubble time for near-circular orbits at the 
virial radius where the tidal field is comparatively weak and 
the mass loss rate is slow. Since the delay depends sensi¬ 
tively on fb, the mass-to-light ratio and the profile of the 
embedded stars we simply note that where it is large is also 
the regime where debris will be streams for tens of Gyr and 
therefore Vstream may be better interpreted as an upper limit 
since some subhalos will not yet have developed visible tidal 
features. Second, since the angles fk e and fk l are propor¬ 
tional to through the energy and angular momentum 

scales their rate of increase is slowed proportionally. The or¬ 
bital width a is not affected and the result > 4^1, still 
holds, therefore shell development is slowed by a factor of 
~ (0.20)"^/® = 1.7. 

4.2.2 Integration limits: approximating observability 

To facilitate a first test of the sensitivity of morphologi¬ 
cal fractions to orbit distributions, we make a number of 
assumptions regarding which mergers will be visible and 
counted as shells or streams, leaving more detailed estimates 
to future work. 

• The threshold /r = 1 is assumed to divide the morpho¬ 
logical classes exactly, independent of orientation. 

• We presume that the existence of coherent debris struc¬ 
tures implies that there has not been a major merger (which 
we take to be mergers greater than 10:1) during disruption 
and therefore and set ^max = 0.1. 

• Host mass growth is neglected, consistent with the pre¬ 
vious assumption. 

• The surface brightness and therefore detectability of 
substructures are a function of mass, orbit and time (John¬ 
ston et al. 2008; Gooper et al. 2010). As an estimate we 
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Zinfall — 0.2115 Zinfnll — 0.5 Zhifall — 1.0 



Figure 9. Slices through the four-dimensional space that Equation 9 integrates over. The outer black line indicates the host halo’s 
escape velocity while the twelve interior lines show where fi = 1, i.e. the transition from shells to streams, in different merger scenarios. 
Any satellite whose orbit places it below the relevant line at the given redshift, host mass, and mass ratio will create debris with a shell 
morphology. The dependence on host mass at fixed ^ is due to the mass-concentration relation changing the slope of the host potential 
at perigalacticon and therefore affecting e^. 


assume that debris structures from 2max = 0.5 (~ 5 Gyr) 
will still have high enough surface brightness to be visible if 
the progenitor was sufficiently massive, choosing sufficiently 
to be ^min = 0.05 for convenient comparison with measured 
infall distributions. Infalls that occur at higher redshift or 
with smaller mass ratios are not counted. A simple extension 
would be a mass-ratio-dependent maximum infall redshift; 
a full accounting would require more complex analysis to as¬ 
sess the influence of debris morphology and orientation on 
the lifetime of observability. 

• The minimum infall redshift for which debris may be 
extended enough to be noticeable is a function of orbital pe¬ 
riod and depends also in fb through the time delay discussed 
in Section 4.2 but for simplicity we choose Zmin = 0.1 (~ 1.3 
Gyr). Another way to define this constraint would be to use 
2min = 0 but Only count mergers with max('If) 5 , \&_l) > 0 
where 0 is some minimum angle; for example, a few times 
the angle subtended by rtide at apogalacticon. 

• Motivated by the results of Section 4.2 we compute the 
scales Ba and la using 20% of the initial satellite mass to cap¬ 
ture the effect of stripping the extended dark matter before 
any stars are removed. 

• Finally, we assume that dynamical friction and as- 
phericity negligibly alter the orbital parameters during the 
erosion of the extended dark halo. 

All of these effects are worthy of further investigation. 

4.3 Sensitivity to orbital infall distributions 

In this section we explore how much influence infall distribu¬ 
tions have on observed number of galaxies that contain shells 
and streams. First we will consider two arbitrarily chosen 
distributions P(Vr, Ve) that are widely separated and have 
simple functional forms (Section 4.3.1). Having established 
this baseline, we consider the differences that result from 
two similar infall distributions measured in the same sim¬ 
ulation with the same technique but at different redshifts 


(Section 4.3.2) and finally to mass dependence in the distri¬ 
butions (Section 4.3.3) using the indicator Naheu{M). 

Evaluating the four-dimensional integral in Equation 9 
numerically is difficult because computing the derivatives of 
orbital quantities that are used to calculate /i requires four 
optimizations and seven integrations. Applying an adaptive 
quadrature method in four dimensions would require tens 
of millions of evaluations and is therefore impractical, so in¬ 
stead we use a significantly modified and parallelized version 
of MCiNT^, a Monte Garlo integrator, to reach acceptable 
~ 1% accuracy with a few tens of thousands of evaluations. 

One way to visualize the division of orbits into shell¬ 
forming and stream-forming is to consider the (Vr, Ve) plane 
and draw lines of /r = 1 (the selected value for the morpho¬ 
logical transition) for different mass ratios and host masses 
at a single redshift. This is shown in Fig. 9. The infall dis¬ 
tribution is then a probability density function (pdf) of the 
coordinate axes and the fractional ratio of shells and streams 
at the given (M, z) is the integral of that pdf on either side 
of the p = 1 line. 


4-3.1 A simple test - distinguishing infall distributions 

As a proof of concept, we use Equation 9 to compute the 
expected fraction of galaxies that would be observed with 
shell or stream debris under the assumptions of Section 4.2.2 
given the two simple, Gaussian infall distributions shown in 
(Vr, Vs) space in the upper panel of Fig. 10. The first, dis¬ 
played with red dashed contours, is axis-aligned and cen¬ 
tered on (Vr,Ve) = (0.7,0.4); this set of orbits is radially 
biased. The second distribution, shown in solid blue con¬ 
tours, is centered on (VrjVs) = (0.75,0.75) and oriented 
at a 45° angle to the axes which produces orbits that have 


^ https://pypi.python.org/pypi/mcint/ 
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Figure 10. Computed ratios of shells and streams given two sim¬ 
ple Gaussian infall distributions. Top: contours containing 5%, 
33%, and 67% of orbits for the selected infall distributions. Dis¬ 
tribution 1 is significantly more radially biased than Distribution 
2 and has lower average energy. Center: Computed fraction of 
galaxies that would be observed to host shells or streams given 
the above distributions. As expected, the more radial orbits of 
Distribution 1 produce more shells. Bottom: the ratio of the num¬ 
ber of shell-hosting to stream-hosting galaxies varies dramatically 
with the avaliable orbital parameters, from approximately 3:4 for 
Distribution 1 to nearly 1:10 for Distribution 2. 


higher average energy and a smaller dispersion in energy 
than the first distribution. 

In this scenario the classical intuition that a more ra¬ 
dial population of satellite orbits will produce more shells 
clearly holds true, as shown in the center panel of Fig. 10. 
The ratio of the expected number of galaxies hosting shells 
to those with streams for the relatively eccentric satellites 
of Distribution 1 is Afgheii/A^stream = 75% for a host mass 
of IO^^Mq, decreasing to 67% for a host mass of 
shells contribute significantly to the global population of de¬ 
bris structures. Conversely, when the satellites are given a 
more equitable distribution of radial and tangential veloci¬ 
ties as in Distribution 2, streams dominate the substructure 
population by nearly 10 to 1. 

This experiment give credence to the idea that the de¬ 
bris structure population contains extractable information 
about the way satellites are accreted and enabling the ex¬ 
clusion large classes of infall distributions. In the following 
sections we consider distributions as measured from within 
a single cosmological model; while the ability to distinguish 
between such similar distributions is more speculative, these 
examples provide a test of the precision necessary from both 
observations and simulations to provide more detailed con¬ 
straints on orbit distributions. 


4-3.2 Comparing infall distributions as measured in 
cosmological simulations 

In Fig. 11 we compare the number of debris structures of 
each type expected from the Benson (2005) distributions as 
measured at z = 0 and z = 1 (i.e. BO and Bl), which are 
shown in the left and center panels, respectively. Both have 
the form of a two-dimensional Maxwell-Boltzmann distri¬ 
bution in tangential velocity with a Gaussian distribution 
of radial velocities but the mean and dispersion in the 14 
Gaussian is a function of Ve to accommodate covariance be¬ 
tween the velocities. Our experiment asks how different the 
number of structures at z = 0 would be if the infall distri¬ 
bution had frozen as it was at z = 1 instead of evolving into 
BO. For Bl the covariance between Vr and Ve is somewhat 
weaker and the distribution overall is more sharply peaked 
at lower values for both velocities. 

The right panel of Fig. 11 shows the computed 
values of Nshel\{M) and Vatream(A7) for Mhoat/MQ £ 
{10^^, 10^®,10^“^}. Because Bl is more strongly peaked, a 
larger fraction of orbits are bound and therefore the total 
number of debris structures is 11% larger. While in both 
cases the structures are dominated by streams, the lower 
mean tangential velocity of Bl produces many more shells; 
averaging over the mass bins, there are 38% more shells but 
only 7% more streams in this case. It is encouraging that 
quite subtle differences in infall distributions actually mea¬ 
sured in simulations result in different shell fractions relative 
to each other and to the number of streams that are expected 
to develop simultaneously. 

4-3.3 Mass-dependent infall 

As noted in Section 4.1.3 above, some authors have found 
that P(Vr,Ve) varies as a function of host halo mass. In 
particular, Jiang et al. (2015) fit the infall distribution in 
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Figure 11. Left: visualization of the BO infall parameter probability density. The black line is the host escape velocity while the color 
and contours encode the probability of an infall at a particular point in this space. Center: same but for Bl. Note the lower average total 
velocity and the sharper peak. Right: Number of shells and streams as computed from Equation 9 using the limits described in the text. 
While streams (triangle markers) dominate, the number of shells (circle markers) increases 38% when changing the infall distribution 
from BO (blue dashed lines) to Bl (red solid lines). 


a 3 X 3 grid of host masses and mass ratios using a Voigt 
profile in total velocity and an exponential in radial velocity; 
we transform to the (Vr,Ve) plane and show the results for 
the highest mass ratio bin, 0.05 < C < 0.5, in Fig. 12. 

To check whether mass dependence will alter A^sheii and 
A^stream we compute them using two scenarios: first assum¬ 
ing that all host haloes, regardless of mass, receive satellite 
infalls as described by the lO^^M©, high-mass-ratio distribu¬ 
tion (J12, leftmost) and secondly assigning them their actual 
host-mass-dependent, high-mass-ratio orbital infall distribu¬ 
tion (JC, the three distributions shown). 

Since the largest mass ratio mergers are more radially 
biased (compare also Tormen 1997; Vitvitska et al. 2002) the 
fractional contribution of shells to the debris population is 
much larger for both J12 and JC (Afgheii/A^stream ~ 40 — 60%) 
than for BO or Bl where all values of ^ are considered 
together, which instead results in Afaheii/Vatream of only 
~ 15 — 20%. Between J12 and JC, the IO^^Mq values are 
identical by construction, however the differences between 
the IO^^Mq and lO^^M© infall distributions increase the 
number of shells calculated around the larger host by 29 
per cent while negligibly altering the number of streams. 


5 DISCUSSION: OBSERVATIONAL 
PROSPECTS 

The estimates of the frequency of shell-like and stream-like 
debris signatures in Section 4 suggest that a survey that 
can accurately measure the prevalence of streams and shells 
would be capable of providing constraints on orbital infall 
distributions. There is one observational sample that is in¬ 
teresting to compare with our preliminary estimates. Atkin¬ 
son, Abraham & Ferguson (2013) surveyed 1700 galaxies to 
27.7 mag/arcsec^ in the g' band, searching for evidence of 
tidal debris. Approximately 1 in 6 galaxies had strong indi¬ 
cations of some type of tidal feature allowing percent-level 
constraints on the fraction of galaxies containing debris of 


several morphological types, as sorted by visual classifica¬ 
tion. 

We can make a preliminary comparison between their 
results and our predictions based on the Jiang et al. (2015) 
infall distributions (Fig. 12). We simplify their six morpho¬ 
logical categories by combining ‘shells’ and ‘fans’ into a sin¬ 
gle population along with ‘streams’ and ‘linear’ features, 
to better match our dichotomy. ‘Diffuse’ and ‘arm’ struc¬ 
tures are presumed to be highly mixed remnants and parts 
of the host galaxy, respectively, and therefore we exclude 
them from our accounting to match the analysis above. By 
assuming that the observed galaxies follow the stellar-to- 
halo mass relation found by Mandelbaum et al. (2006) and 
combining the observed stellar masses into 0.4 dex bins, the 
Atkinson, Abraham & Ferguson (2013) results suggest that 
at host halo masses of IO^^Mq approximately 3 per cent 
of the observed galaxies have shells and 6 per cent have 
streams, while of those nearer IO^^Mq 8 per cent have shells 
and 10 per cent have streams. Comparing with Fig. 12, one 
sees that the analytic estimates overpredict both types of 
debris structures by a factor of ~ 2. In addition, the ratio of 
shells to streams is higher in the observations and is more 
strongly dependent on mass. While it may be tempting to in¬ 
terpret these numbers as indications that the model merger 
rate is too high with orbits that are excessively tangentially 
biased, several questions remain. These differences could in¬ 
stead be attributed to the simplifications made in the inte¬ 
gration limits that were chosen to mimic a surface brightness 
limits in the observations or in differences between their vi¬ 
sual classification of debris types and our analytic criteria 
for splitting shells and streams. More carefully matching the 
surface brightness cut to the survey parameters will allow us 
to probe these questions. 

It is clear that extensive work is necessary before a 
conclusive comparison between observational and simula¬ 
tion can be made. This requires (i) an unbiased, uniform 
survey containing at least a few hundred galaxies with de¬ 
tected debris structure, such as the Atkinson, Abraham 
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Figure 12. Infall distributions and resultant average number of shells and streams from J12 and JC. Left panels: The transformed 
velocity distributions in the highest mass ratio bin from Jiang et al. (2015). There is a trend of decreasing covariance between radial 
and tangential velocity with increasing host mass as well as lower average tangential velocity at larger host masses. Right: the resultant 
debris populations, first when assuming all hosts receive infalls according to the lowest host mass infall distribution (J12, gray lines) 
and then assigning each their measured probabilities (JC, black lines). The same values are found In each case at Muost = IO^^Mq by 
construction. 


& Ferguson (2013) sample; (ii) a set of satellite disrup¬ 
tion simulations with different mass ratios and orbits type, 
realistically-embedded stellar distributions, considering all 
viewing orientations and including contamination from the 
host galaxy’s light; and (hi) an objective and automated 
method for classifying features to a given surface brightness 
limit in both simulations and observations in the same way. 

Looking to the future, the Large Synoptic Survey Tele¬ 
scope will probe a region ~ 20000 deg^ and ultimately reach 
a surface brightness sensitivity of r ~ 29 mag arcsec”^ 
(Ivezic et al. 2008). Such a deep, uniform survey will un¬ 
veil vast numbers of substructures. This prospect provides 
strong motivation for addressing the second and third re¬ 
quirements so that the orbits of infalling structures can be 
recovered. 


6 CONCLUSION 

In this work we have studied the formation of the major mor¬ 
phological classes of tidal debris from minor galactic merg¬ 
ers - shells and streams - and found that a simple physical 
model of stream production based on offsets in energy be¬ 
tween the unbound particles and the progenitor satellite is 
extensible to shells if the consequences of angular momen¬ 
tum offsets are also included. This is quantified through a 
morphology metric which is shown to provide a good match 
to visual classification. For an individual merger the result¬ 
ing morphological class is intrinsically time-dependent be¬ 
cause the energy effects are bounded by the size of a sin¬ 
gle orbit but those of the angular momentum are not. We 
demonstrate how the distribution of merging satellite or¬ 
bital parameters measured in cosmological N-body simula¬ 
tions can be interpreted through a morphology metric as 
predictions for the fraction of galaxies that have experienced 
stream- or shell-creating mergers. Our results show that dif¬ 
ferent infall distributions produce results that are plausibly 
distinguishable using studies of low-surface-brightness fea¬ 
tures around nearby galaxies. Tidal debris morphology thus 
provides unique access to orbit distributions - an as yet un¬ 
explored part of the accretion history. 
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